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Although the rutile structure of Ti02 is stable at high temperatures, the conventional quasiharmonic 
approximation predicts that several acoustic phonons decrease anomalously to zero frequency with 
thermal expansion, incorrectly predicting a structural collapse at temperatures well below 1000 K. 

Inelastic neutron scattering was used to measure the temperature dependence of the phonon density 
of states (DOS) of rutile Ti02 from 300 to 1373 K. Surprisingly, these anomalous acoustic phonons were 
found to increase in frequency with temperature. First-principles calculations showed that with lattice 
expansion, the potentials for the anomalous acoustic phonons transform from quadratic to quartic, 
stabilizing the rutile phase at high temperatures. In these modes, the vibrational displacements of 
adjacent Ti and O atoms cause variations in hybridization of 3d electrons of Ti and 2p electrons of O 
atoms. With thermal expansion, the energy variation in this "phonon-tracked hybridization" flattens 
the bottom of the interatomic potential well between Ti and O atoms, and induces a quarticity in the 
phonon potential. 


I. INTRODUCTION 

Titanium dioxide (Ti02) is of longstanding interest in 
physics, chemistry, surface science and materials sci¬ 
ence. It has numerous technological applications in¬ 
cluding pigments, optical coatings, catalysis, solar cells 
and gas sensorsii’^ Rutile is the stable phase of Ti02 at 
temperatures below 1800 K, and the two other naturally- 
occurring phases of Ti 02 , anatase and brookite, both con¬ 
vert to rutile upon heating;^ The thermal properties of 
rutile Ti 02 are central for applications that involve heat 
generation, thermal transport and temperature-driven 
phase transitions. Nevertheless, even after numerous 
experimental^— and theoreticalii“— investigations, the 
stability of the rutile phase at high temperatures remains 
poorly understood. 

The quasiharmonic approximation (QHA) is based on 
how phonon frequencies change with volume. In the 
QHA, all shifts of phonon frequencies from their low 
temperature values are the result of thermal expansion 
alone— i2i For a lattice expansion of 0.5%, the energies 
of transverse phonons decrease by 10% to 50% in the 
QHA, giving Griineisen parameters as large as 100 for 
rutile TiQ 2 — ^ In the QHA, many transverse acoustic 
modes decrease to zero frequency at a lattice expansion 
of 0.75%, so the QHA predicts a collapse of the rutile 
structure at temperatures below 1000 K. The QHA does 
not account for all non-harmonic effects, however. Al¬ 
though the QHA accounts for some frequency shifts, 
the phonon modes are still assumed to be harmonic, 
non-interacting, and their energies depend only on the 
volume of the crystal. 

In non-harmonic potentials, phonon-phonon interac¬ 
tions are responsible for pure anharmonicity that short¬ 
ens phonon lifetimes and shifts phonon frequencies. 
Anharmonicity competes with quasiharmonicity to al¬ 


ter the stability of phases at high temperatures, as has 
been shown, for example, with experiments and frozen 
phonon calculations on bcc Zr— and the possible stabi¬ 
lization of bcc Fe-Ni alloys at conditions of the Earth's 
core— For PbTe and ScFa, there are recent reports of 
anharmonicity being so large that both the QHA and 
anharmonic perturbation theory fail dramaticallyi^i^^ 
These cases are suitable for ah initio molecular dy¬ 
namics (AIMD) simulations, however, which naturally 
includes all orders of anharmonicityi^— The AIMD 
method should be reliable when the electrons are near 
their ground states and the nuclear motions are classical. 

In this work, we study phonons at high temperatures 
in rutile TiQ 2 by neutron inelastic scattering, AIMD, and 
other methods. We report a remarkable quartic anhar¬ 
monicity that develops because the hybridization be¬ 
tween electron states at Ti and Q atoms changes dy¬ 
namically during atom vibrations, and these changes 
are highly sensitive to lattice expansion. This "phonon- 
tracked hybridization" induces the phonon quarticity 
that is essential for stabilizing the rutile phase at high 
temperatures, and affects properties such as ferroelec- 
tricity and thermal transport. 


II. EXPERIMENT 
A. Inelastic neutron scattering 

Samples were commercial TiQ 2 powder (Alfa Ae- 
sar. Ward Hill, MA) with a rutile phase fraction of at 
least 99.9%. The sample powder was loaded into a 
niobium sachet for neutron scattering measurements 
at high temperatures. Measurements were performed 
with the ARCS spectrometer at the Qak Ridge National 
Laborator}*^^, using incident neutron energies of 75 meV 
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FIG. 1. S(Q, E) spectra of rutile Ti02 measured at (a) 300K and (b) 1373K, respectively. The incident energy of neutron beam is 30 
meV. 


and 30 meV. A vacuum furnace with niobium radiation 
shields was used for measurements at 300, 673, 1073, 
and 1373 K. Furnace backgrounds were measured with 
empty sample containers at each temperature. Data 
reduction was performed with the software package 
DRCS3 The raw data of individual neutron detection 
events were first binned into scattering angle and energy 
transfer, E, and normalized by the proton current on tar¬ 
get. The data were corrected for furnace background and 
detector efficiency. The data were then rebinned into in¬ 
tensity, S(Q, E), where fiQ is the momentum transfer to 
the sample. This S(Q, E) is averaged over all crystal¬ 
lographic directions, and is a reasonably well-balanced 
sampling of all phonons in the material. 

B. Results 

Surprisingly, temperature causes a global increase in 
energy of the entire TA branch. The lowest phonon dis¬ 
persions formed by the TA branch at 14 meV at 300 K 
increase significantly in energy at 1373 K, as seen in the 
S(Q, E) phonon spectra (Fig. U). There are no disper¬ 
sions that soften between 2 and 4 A“^, and there are 
no signs of phonons that collapse to zero frequency. 
This anomalous thermal stiffening is also apparent in 
the phonon DOS curves, which were obtained after cor¬ 
rections for multiphonon and multiple scattering using 
the getdos packagei^ As shown in Fig. |2la), peak 1 (cen¬ 
tered at 14 meV) exhibits an unusual stiffening as large 
as 2.7meV, and it remains sharp. The large stiffening of 
this DOS peak is even more apparent in measurements 
with higher energy resolution, using an incident energy 
of 30 meV (see Fig.|3la)). On the other hand, as shown in 
Fig.lUa) and Fig.|3la), most other features of the phonon 
DOS above 20meV undergo substantial softening and 
broadening with temperature, in good agreement with 
the individual phonons measured by Raman spectrom¬ 
etry at high temperatures It is well-known that rutile 
Ti02 has a ferroelectric soft mode A 2 ,, at the F-point that 
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FIG. 2. (a) Neutron weighted phonon DOS of rutile Ti02 from 
measurements at temperatures from 300 to 1373 K with an in¬ 
cident energy of 75 meV. Five vertical lines are aligned to peak 
centers at 300 K. (b) Phonon DOS of rutile Ti02 at temperatures 
of 300 and 1373 K calculated with first principles MD simula¬ 
tions and the FTVAC method. 
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FIG. 3. (a) Neutron weighted phonon DOS of rutile Ti02 from 
measurements at temperatures from 300 to 1373 K with an inci¬ 
dent energy of 30 meV (so data above 26 meV are less reliable). 
The dashed spectrum corresponds to the experimental result 
at 300 K, shifted vertically for comparison at each temperature, 
(b) Simulated total phonon DOS from the QHA at 300 K (black) 
and 1373 K (red). 


stiffens with heatingi^ but the global stiffening of the 
entire lower-lying TA branch is surprising. 

Figure|3[b) presents the DOS calculated with the quasi¬ 
harmonic model at 1373 K, showing that contrary to the 
experimental trend in Fig. |3la), the first peak in the 
phonon DOS undergoes a large softening with temper¬ 
ature. Consistent with previous results/i^^ our QFIA 
calculations predict that modes in the TA branch soften 
through zero frequency with the thermal expansion of 
1373 K, giving imaginary frequencies in the DOS (shown 
as negative in Fig. |3|b)) that would destabilize the rutile 
structure at high temperatures. 


TABLE I. Lattice parameters and linear expansivities (10 ^/K) 
measured by experiment and calculated by MD simulation 



a (A) 

(Xa 

c(A) 

ac 

U 

Exp(300K) 

4.592 


2.958 


0.3092 

Exp(1373K) 

4.635 

8.73 

2.993 

11.0 

0.3123 

Cal(300K) 

4.565 


2.935 


0.3058 

Cal(1373K) 

4.603 

8.75 

2.967 

10.2 

0.3087 


III. CALCULATIONS 

A. First principles molecular dynamics simulations and 
techniques for calculating the vibrational spectra 

First-principles calculations using the local density ap¬ 
proximation (LDA) of density functional theory (DFT)^ 
were performed with the VASP package<^2i2i Projector 
augmented wave pseudopotentials and a plane wave 
basis set with an energy cutoff of 500 eV were used in 
all calculations. Previous work showed that for best ac¬ 
curacy, the Ti pseudopotential should treat the semicore 
3s and 3p states as valence /^'^^d^ and we took this ap¬ 
proach with a similar LDA functional. Our calculated 
elastic properties, lattice d 5 mamics, and dielectric prop¬ 
erties derived from the optimized structure for 0 K, were 
in good agreement with results from experiment and 
from previous DFT calculations. Phonon dispersions 
and phonon DOS spectra were obtained in the quasihar¬ 
monic approximation by minimizing the vibrational free 
energy as a function of volume 

+00 

F{V, T) = Eo + Jg(aj) + hTln(l - J do; ,(1) 

— CO 

where Eq is the energy calculated for the relaxed struc¬ 
ture at T = 0 K. 

First-principles Born-Oppenheimer AIMD simula¬ 
tions for a 2 X 2 X 4 supercell and a 2 X 2 x 1 fc-point 
sampling were performed to thermally excite phonons 
to the target temperatures of 300 and 1373 K. For each 
temperature, the system was first equilibrated for 3 ps as 
an NVT ensemble with temperature control by a Nose 
thermostat, then simulated as an NVE ensemble for 20 ps 
with time steps of 1 fs. Good relaxations with residual 
pressures below 0.5 GPa were achieved in each calcu¬ 
lation that accounted for thermal expansion. The lattice 
parameters and thermal expansivities determined by the 
simulations are presented in Table |T] 

1. Fourier transformed velocity autocorrelation method 

Velocity trajectories were extracted from the MD sim¬ 
ulation at each temperature, and were then transformed 
to the corresponding vibrational energy and/or momen¬ 
tum domain ] i9r^6-28,35 Because the FTVAG method does 
not assume a form for the Flamiltonian, it is a robust tool 
for obtaining vibrational spectra from MD simulations, 
even with strong anharmonicity. The phonon DOS is 
given by 

y(co) = 2], re“'“V„,b(f)4o(0)>df, (2) 

n,h 

where () is an ensemble average, and is the ve¬ 

locity of the atom b in the unit cell n at time t. Further 
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FIG. 4. Diffuse curves are TDEP phonon dispersions below 25 meV at (a) 300 K, and (b) 1373 K, compared with the results from 
the FTVAC method (red circles). The white curves are phonon dispersions for the quasiharmonicity plus quartic anharmonicity 
calculated with all set to zero in Eq.|3 In (b), the dispersions are compared to the quasiharmonic dispersions (orange dashed 
curve) and the single quartic oscillator model (orange triangles), (c) Frozen phonon potential (black) of TA mode at R point with 
q = (0.5,0,0.5) at 1373 K, showing the harmonic component (red) and quartic component (blue). The low temperature potential 
surface is also shown (dashed black). 


TABLE II. Frequencies (meV) of the transverse acoustic modes of rutile at ambient and high temperature, from experiment and 
from calculations with QHA, MD and TDEP methods 



X 

(0.5,0,0.25) 

R 

Z 

(0,0,0.25) 

F 

M 

(0.5,0.5,0.25) 

A 

T = 300K 

Exp^ 

12.77 



15.25 


0 

12.15 



QHA^ 

11.04 

17.30 

10.42 

10.42 

7.15 

0 

9.55 

5.65 

9.92 

QHA^ 

13.52 


13.02 

12.90 


0 

9.05 


12.03 

QHA 

11.68 

17.67 

12.30 

12.60 

8.33 

0 

9.15 

7.09 

12.09 

MD 

13.60 

19.88 

13.49 

15.36 

10.36 

0 

10.95 

9.90 

13.07 

TDEP 

13.60 

22.27 

8.51 

16.26 

10.40 

0 

11.01 

6.01 

9.20 

T = 1373K 

QHA^ 

5.81 

12.85 

4.09 

7.61 

3.06 i 

(>7)i 

4.28 

5.75 i 

8.63 

QHA 

3.13 

12.59 

2.02 

9.07 

3.78 i 

9.83 i 

3.50 

5.82 i 

9.83 

MD 

15.41 

19.81 

16.93 

18.04 

13.34 

0 

12.27 

13.78 

14.28 

TDEP 

16.91 

21.12 

13.53 

19.63 

13.35 

0 

11.07 

11.48 

11.57 


® Experimental data from Ref. @| 

’’ First principles calculations from Ref. jl^ 

' First principles calculations from Ref. Il6ll 

First principles calculations from Ref. with ambient temperature but 1% isotropic lattice expansion 


projection of the phonon modes onto each k point in the 
Brillouin zone was performed by computing the phonon 
power spectrum with the FTVAC method, with a res¬ 
olution determined by the size of the supercell in the 
simulation. 


2. Temperature-dependent effective potential method 

In general, the cubic phonon anharmonicity con¬ 
tributes to both the phonon energy shift and the life¬ 
time broadening, whereas the quartic anharmonicity 


contributes only to the phonon energy shifti^^^ To dis¬ 
tinguish the roles of cubic and quartic anharmonicity, 
the temperature-dependent effective potential (TDEP) 
method^Si^S was used. In the TDEP method, an effec¬ 
tive model Hamiltonian is used to sample the potential 
energy surface, not at the equilibrium positions of atoms, 
but at the most probable positions for a given tempera- 
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ture in an MD simulations^ 

i ija^ 

+1 E ' 


a p 
U. 
1 


ijkafiy 


( 3 ) 


where cpij and are second- and third-order force con¬ 
stants, p is momentum, and u“ is the Cartesian compo¬ 
nent a of fhe displacemenf of atom i. In fhe fiffing, fhe 
"effective" harmonic force consfanfs cpij are renormalized 
by fhe quarfic anharmonicify. The cubic anharmonic- 
ify, however, is largely accounted for by fhe fhird-order 
force consfanfs ^pip and can be undersfood in terms of 
fhe fhird-order phonon self-energy fhaf causes linewidfh 

broadeningiSSdS 

The resulfing Hamilfonian was used to obfain fhe 
renormalized phonon dispersions (TDEP specfra) ac- 
counfing for bofh fhe anharmonic shiffs A, and broad- 
enings T, of fhe mode qj. These are derived from fhe 
real and imaginary parfs of fhe cubic self-energies 
respecfivelyjS^ 


A{qj-, ^) = -^ E E +q 2 -^ 

qih Hih 

r ni+n 2 + l n\ + n 2 + l 

+ (Ui -(- Cj02 — CjOi — (i)2 

»1 - »2 _ Wi - W2 1 

Q — <x)i -t (x>2 + <jJl ~ OJ 2 

r(^y;C^) = ^ E E A(ifi + q 2 -^ 

qih qih 

x[(ni -h n 2 + 1) 6(C2 - coi - 0 ) 2 ) 

+ 2(ni — M 2 ) 6(Q -H coi — CiJ2)j > (5) 

where Q is fhe renormalized phonon frequency and p de- 
nofes fhe Cauchy principal parf. The V(.)'s are elemenfs 
of fhe Fourier fransformed fhird order force consfanfs 
pJijk obfained in fhe TDEP mefhod. The A{qi + q 2 — ^ 
ensures conservafion of momenfum. 


fhe eigenvectors of fhe mode. Self-consisfenf wavefunc- 
fions were fhen computed and projecfed onfo a local 
basis of Slafer t 5 rpe orbifals wifh proper orfhonormaliza- 
fion, giving fhe COEIP specfra of fhe system. Since nega- 
five values of fhe COEIP suggesf bonding confribufions, 
- COHP specfra are presented, following convenfion. 


C. Results 

As shown in Fig. |2|b), the neutron-weighted DOS 
spectra calculated with FTVAC method reproduced the 
thermal shifts and broadenings observed in the exper¬ 
imental DOS of Fig. Ilia). In particular, upon heating 
fo 1373 K, an anomalous stiffening of peak 1 by 2.1 meV 
was obfained. 

Figure |4] shows fhe vibrafional energies of fhe TA 
branch, calculated by fhe FTVAC mefhod wifh AIMD 
frajecfories. From 300 fo 1373 K, fhe TA branch increases 
in energy by an average of abouf 2.1 meV. Especially for 
fhis TA branch. Fig. lUb) shows an enormous discrep¬ 
ancy of phonon energies befween fhe MD calculafion 
and fhe QHA (orange dashed line) af 1373 K. The unsfa- 
ble phonon modes predicted by fhe QFIA are fully sfable 
in fhe AIMD simulations af high femperafures, however. 

Using fhe same MD frajecfories as for fhe FTVAC 
mefhod, fhe calculated TDEP dispersions agree well 
wifh fhe FTVAC resulfs as shown in Fig. Illa)(b), and 
also in Table HIl (The full TDEP dispersions are shown in 
fhe Supplemenfal Maferial4^) The cubic anharmonicify 
of rutile Ti 02 is sfrong for phonons af energies above 
25 meV/i^ causing broadening of fhe phonon DOS of Fig. 
|2]and broadening of fhe dispersions^^ Neverfheless, af 
1373 K fhe TA modes below 20meV have only small 
linewidfh broadenings. Furfhermore, fhey are close in 
energy fo fhose calculated if all ipijk are sef fo zero in Eq. 
13 showing fhe dominance of quarfic anharmonicify and 
fhe small cubic anharmonicify of fhe TA modes. 


IV. DISCUSSION 
A. Phonon quarticity of TA branch 


B. Crystal orbital Hamilton population method 

The Crystal Orbital Flamilton Population Method 
(COEIP) projects plane waves to local orbital basis func¬ 
tions and is used to partition the band-structure energy 
into the bonding, nonbonding and antibonding interac¬ 
tions between two atoms. COHP provides a semiquan- 
titative interpretation of fhe nef bonding characferisfics. 

In fhis work, COHP specfra were calculated wifh fhe 
package Focal-Orbifal Basis Suife Towards Elecfronic- 
Sfrucfure Reconsfrucfion (FOBSTER)4id^ The phonon 
displacemenf patterns were generated in a 2 x 2 x 4 su¬ 
percell wifh a 8 X 8 X 4 fc-poinf sampling according fo 


For more defails abouf fhe anomalous anharmonicify 
of fhe TA modes, fhe frozen phonon mefhod was used fo 
calculafe pofenfial energy surfaces for specific phonons, 
as in Fig. SJc). Af 300 K fhe pofenfial energy of fhe TA 
mode af fhe R-poinf is nearly quadratic, wifh a small 
quarfic parf. Wifh fhe laffice expansion characferisfic of 
1373 K, fhe pofenfial energy curve fransforms fo being 
nearly quarfic. 

In facf, for all modes in fhe TA branch fhaf were eval- 
uafed by fhe frozen phonon mefhod, fhe pofenfial en¬ 
ergy surface develops a quarfic form wifh laffice expan¬ 
sion (see Supplemenfal Material!^). Even more inferesf- 
ing behavior occurs along fhe directions of M-A, T-Z 
and T-X, where fhe pofenfial energy surface develops a 
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FIG. 5. (a) Displacements of atoms for the TA mode at the R-point in the (1-10) plane. Light blue spheres are Ti atoms, and O 
atoms are red. Arrows depict distortions of the structural units (dashed rhombuses). The rotational movements of structures, or 
the "ring" patterns, are indicated with circled arrows, (b) ELF isosurface of a "ring" shown in (a) with the value of 0.3. The ELF 
increase is apparent in the bond of shorter distance owing to the ring displacement. The ELF is the probability measure of finding 
an electron at a location given the existence of neighbouring electrons. It ranges from 0 (no electron) to 1 (perfect localization), (c) 
COHP analysis of Ti-O bonds for equilibrium lattice parameter at T=0K (0%), and for linear expansions of 1% and 2%. Shown 
in color are -COHP results for the same structures with the phonon of panel (a) having 0.14 A normal displacements of Ti atoms. 
Curves for 1% and 2% expansion are offset by 0.6 and 1.2. 


Landau-t 5 rpe form, with negative curvature at zero dis¬ 
placement, but a positive quartic shape at large displace¬ 
ments. (The height of the energy barrier between the 
two minima is lower than k^T, however.) For a quantum 
quartic oscillator, the vibrational frequency stiffens with 
temperature owing to the increasing spread between the 
energy levelsi^^^ We assessed a high temperature be¬ 
havior by assigning Boltzmann factors to the different 
oscillator levels derived from frozen phonon potentials, 
giving the energies of the quartic TA modes at 1373 K. 
As shown in Fig. IHb), they are reasonably close to the 
FTVAC and TDEP results. 


B. Phonon-tracked hybridization and its change with 
thermal expansion 

The patterns of atom displacements in the anomalous 
modes at F, R, and along Z-F, F-M, and M-A were iden¬ 
tified, and those for the R point in Fig. [Sja) are t 5 rpical. 
In these anomalous modes, the O atoms were approxi¬ 
mately stationary, and each O atom has a Ti neighbor that 
moves towards it and another Ti neighbor that moves 
away from it by approximately the same amount. These 
modes have "ring" patterns in which displacements of 
Ti atoms rotate a structural unit, and all the O atoms 
see approximately the same change in their Ti neigh¬ 
borhood. In the positive and negative displacements 
of these modes, the O atoms show an accumulation of 
charge in the Ti-O bond of shorter distance and a de¬ 
pletion in the bond of longer distance, as indicated by a 
much higher value of the electron localization function 
(ELF)i^ for the short bond shown in Eig.|3[b). 

We calculated the "bond-weighted" electron DOS by 
partitioning the band structure energy into bonding and 


no-bonding contributions and obtaining the crystal or¬ 
bital Flamilton population (COFIP) spectrum^id^ of ru¬ 
tile with different lattice parameters in a 2 x 2 X 4 super¬ 
cell. Pigurel^c) shows the COHP spectrum of the bonds 
formed by the Ti3d and 02p orbitals between 5.5 and 
3.5 eV below the Permi energy. With lattice expansion (of 
1% or 2%), these bonding states become less favorable, 
and their COHP spread narrows. Also shown in color 
in Pig. |3b) is the COHP with a frozen phonon mode 
at the R-point having 0.14 A displacements of Ti atoms. 
On the scale of thermal energies, the broadening effect 
from the phonon changes considerably with lattice ex¬ 
pansion. Por the equilibrium lattice parameter (labeled 
"0%"), bonding states at both the top and bottom of the 
Ti 3d - O 2p band are shifted upwards by the phonon. 
With larger lattice parameter, however, the changes take 
place mostly in the states at the bottom of the band, 
which broadens the COHP spectra and helps to lower 
the bonding energy. Additional bonding states also ap¬ 
pear near -3 eV. The overall effect is that the potential 
energy is relatively insensitive to small displacements, 
and the phonon potential has a much flatter bottom with 
1% lattice expansion as shown in Pig.lHc). Por a 2% lin¬ 
ear expansion, a Landau-type phonon potential evolves, 
with local minima away from zero displacement owing 
to the spread of states to the lowest part of the band. 

The effect can be understood by adapting a hard 
sphere model for the chain of -Ti-O-Ti-O- atoms along the 
direction "A" or around the "ring" units in the phonon 
displacement pattern in Pig. ISja). With an increase in 
lattice parameter, the longer Ti-O bond makes a smaller 
contribution to the interatomic force during its vibra¬ 
tional cycle. The shorter Ti-O bond gives a stronger hy¬ 
bridization of Ti 3d - O 2p orbitals as the Ti atom moves 
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FIG. 6. Phonon potentials for the "ring" pattern of Fig. EJa). 
The solid curve is the potential of the R-point mode at 1373 K as 
in Fig.lHc). The dashed curve is the potential of the same mode 
1373 K, but with the "ring" pattern broken by immobilizing the 
Ti atoms along the A-direction in Fig. Ela) . The shape of dashed 
curve is quadratic. 

closer to the O atom. The hybridization serves to offset 
the energy of short-range repulsion. With lattice expan¬ 
sion, the short-range repulsion is weaker, and hybridiza¬ 
tion favors electrons between the shorter Ti-O pairs in 
the phonon displacement pattern. The "ring" patterns 
of the phonons play an important role in increasing the 
degree of hybridization as they complete the electron 
back-donation cycles from the O to the Ti atoms. For ex¬ 
ample, a 3% decrease of Bader effective charge (+2.22 at 
equilibrium) was found for the Ti atoms with 0.14 A dis- 
placments in the "ring" patterns, which is comparable 
to the charge decrease of the Ti atoms during the ferro¬ 
electric transition of rutile However, if Ti atoms along 
the direction "A" were locked down at their equilibrium 
positions so the ring motion is broken, the resulting de¬ 
crease of the effective charge dropped by 50% to 70% in 
the "ring" patterns. The potential was found to rise, and 
was largely quadratic even at 1% or 2% expansion, as 
shown in Fig. |6] 

A macroscopic elastic response to this phonon can also 
be identified with the assistance of Fig. |5la). In equilib¬ 
rium, the apex angles of the rhombuses are all equal, but 
with the rotation by b, the vertical stretching of rhom¬ 
buses along the line A is 2a sin(0 + b), and the contraction 
is 2flsin(6 - b), where 6 is the semi-angle of the rhom¬ 
bus. For small b, a Taylor expansion gives a net vertical 
(or horizontal) distortion of -2ab^ sin 6 (or -2ab^ cos 0). 
The distortions are proportional to b^, while the atom 
displacements in this TA mode are proportional to b. A 


strain energy that goes as the square of this distortion is 
consistent with a quartic potential. 

The hybridization in the Ti-O bond is very sensitive 
to interatomic distance, much as has been noticed in the 
ferroelectric distortion of BaTiOsi^ For rutile Ti02, how¬ 
ever, the hybridization follows the atom displacements 
in thermal phonons (instead of a displacive phase transi¬ 
tion), and this "phonon-tracked hybridization" changes 
with lattice parameter. It provides a source of extreme 
phonon anharmonicity, but also provides thermody¬ 
namic stability for rutile Ti02. It may occur in other 
transition metal oxides that show unusual changes of 
properties with lattice parameter or with structure, and 
such materials may be tunable with composition or pres¬ 
sure to control this effect. Besides altering thermody¬ 
namic phase stability, properties such as ferroelectricity 
and thermal transport will be affected directly. 

V. CONCLUSION 

Several inter-dependent methods proved useful for as¬ 
sessing the anharmonic phonon d 5 mamics of rutile Ti 02 . 
Ab initio methods of molecular d 5 mamics, temperafure- 
dependent effective potential, and frozen phonons gave 
accurate accounts of the thermal stiffening of the TA 
phonons that was measured by inelastic neutron scat¬ 
tering on rutile Ti02. These methods avoided the struc¬ 
tural instability predicted by the quasiharmonic approx¬ 
imation (QHA). With thermal expansion, many acoustic 
modes were unstable in the QHA, as observed previ¬ 
ously. From AIMD simulations and frozen phonon cal¬ 
culations, it was found that when the rutile structure 
undergoes thermal expansion, a hybridization between 
Ti and O atoms gives low-energy electron states when 
the Ti and O atoms are in closest proximity during their 
vibrations. This phonon-tracked hybridization flattens 
the bottom of the potential of the anomalous phonon 
modes, giving quartic potentials that stabilize the rutile 
structure at high temperatures. 
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Supplemental Materials: 

Phonon quarticity induced by changes in phonon-tracked hybridization during lattice 

expansion and its stabilization of rutile TiOi 



FIG. 7. (a) Calculated phonon dispersion of rutile Ti02 along the high symmetry directions, (b) The corresponding phonon DOS. 
The DOS peaks are labeled by the same numbers as those in Fig.l of the paper 



FIG. 8. TDEP phonon spectra along high symmetry points at (a) 300 K and (b) 1373 K 

























Eti*w4tniV) E^Mtav^mVO EnBmrlraoV} 


10 



DIvtauttMnl (71 



aiplflMn)fi< (JJ 



-«.4-0;34.2-a.1 0.D 0.1 0.3 0.3 0.4 
ClItpiKWrwUK?] 



Ot«p4i«vnififi C?> 





FIG. 9. Frozen phonon potentials (black) of TA modes at high symmetry points at 1373 
harmonic component (red) and quartic component (blue). 
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